#include "utils.hpp"
#include "Misc/Laplace.hpp"
#define INF 2 
#define PROP_MAX 0.85 
#define PROP_MIN 0.05

class LogitProp : public Distribution
{
	public:
	LogitProp(){}
	LogitProp(int n)
	{
		_n=n;
		boost::random::uniform_01<> N;
		runif= new RandomG::Random<boost::random::uniform_01<> >(N);
		dn= new boost::math::normal_distribution<>(0,1);
		L=new Laplace;
		_s=mat(n,n);
		_s.fill(1);
		mat B(n,n);
		B.fill(0);
		_B=B;
		
	}
	virtual ~LogitProp()
	{
		delete runif;
		delete L;
		delete dn;
	}

	mat r(int j)
	{
		double sum;
		double u=0;
		mat R(_n,1);
		R.fill(0);
		vector<double> P;
		for(int i=0;i<_n;i++)
		{
			sum=0;
			for(int k=0;k<=i;k++)
			{
				if(i!=k){
					sum+=(_B)(i,k)*R(k,0);
				}else{
					sum+=(_B)(i,k);
				}

			}
			
			double p=cdf((*dn),sum);
			P.push_back(p);		
			u=(*runif)();	
			if(u<p)
			{
				R(i,0)=1;
			}else{
				R(i,0)=0;
			}

		}
		return R;
	}


	//ALWAYS BE SURE THAT r and d a performes sequentially;
	double 	d(mat x, int b=0,int j=0){
		double q=1;
		double p=0;
		double sum;
		x=x.t();
		
		for(int i=0;i<_n;i++)
		{
			sum=0;
			for(int k=0;k<=i;k++)
			{
				if(i!=k){
					sum+=(_B)(i,k)*x(0,k);
				}else{
					sum+=(_B)(i,k);
				}

			}

			p=cdf((*dn),sum);
			if(x(0,i)==1)
			{
				q=q*p;
			}else{
				q=q*(1-p);
			}

		}
		if(b==1)
		{
			return log(q);	
		}else{
			return q;
		}
	}
	mat Var(mat X){
		_X=X;
		int u=0;
		mat X1;
		mat X2;
		mat test=sum(X,0)/X.n_rows;
		uvec _index1= find(test > PROP_MAX);
		uvec _index2= find(test < PROP_MIN);
		uvec _index=merge<uvec>(_index1,_index2);
//		cout << _index1;
//		cout << _index;
		//cout << _index2;
	
		for(int i=0;i<_n;i++)
		{
			//for(int k=0;k<i;k++)
		//	{
			mat X2=X(span::all,i);
			double prop=as_scalar(sum(X2))/X.n_rows;
			if( (prop>PROP_MIN) & (prop<PROP_MAX))
			{	
				mat D(X.n_rows,i+1);
				mat B(1,i+1);
				if(i!=0)
				{
					mat X1=X(span::all,span(0,i-1));
					mat B1=_B(i,span(0,i));
					u=0;
					for(int jj=0;jj<_index.n_elem;jj++)
					{
						if(_index(jj)<=i-1)
						{
							mat temp=del_one(X1,_index(jj)-u);
							X1=temp;
							mat temp2=del_one(B1,_index(jj)-u);
							B1=temp2;
							u++;
						}					
					}
					if(X1.n_cols<1)
					{
						D.resize(X.n_rows,1);
						D.fill(1);
						B.resize(1,1);
						B.fill(_B(0,0));
					}else{
						B=B1;
						D=add_one(X1,1);	
					}
				}else{
					D.fill(1);
					B.fill(_B(0,0));
				}
				double *Y=new double[X2.n_rows];
				mat2vector<double>(X2,Y);
				//cout << D;
				//for(int g=0;g<X.n_rows;g++){cout << Y[g];}
				try{
				 	(*L)(D,Y,_m,(_s)(0,0),B.t());
				}
				catch(...)
				{
					cout <<  "\nModel under consideration for proposal is singular"; 
				}
				mat Foo=((*L).Get_Mu()).t();
				//cout << "Foo " << Foo;
				int g=0;
				for(int k=0;k<=i;k++)
				{
					int foo=is_present<uvec>(_index,k);
					int foo2=is_present<uvec>(_index2,k);
					if((foo==0) & (foo2==0))
					{
						(_B)(i,k)=Foo(0,g);
						g++;
					}else{
						(_B)(i,k)=0;
					}

				}
				//cout << "_B " << _B;
				delete[] Y;
			}else if(PROP_MIN<prop){
				for(int k=0;k<=i;k++)
				{
					if(k==i)
					{
						if(prop!=1)
						{
							(_B)(i,k)=quantile(*dn,prop);
						}else{
							(_B)(i,k)=INF;
						}
					}else{
						(_B)(i,k)=0;
					}

				}
			}else{
				for(int k=0;k<=i;k++)
				{
					if(k==i)
					{
						if(prop!=0)
						{
							(_B)(i,k)=quantile(*dn,prop);
						}else{
							(_B)(i,k)=-INF;

						}
					}else{
						(_B)(i,k)=0;
					}

				}
			}
		//	}
		}
	//	cout << _B;
		return _B;
		//gets B
	}
	mat scrambled(int N)
	{
		mat X(N,_n);
		for(int i=0;i<N;i++)
		{
			X(i,span::all)=this->r(i);
		}
		return X;
	}
	void Set_n(int n){		}
	void Set_m(mat m){		}
	mat Get_m(void){return _m;}
	mat Get_s(void){return _s;}
	void Set_s(mat s){	_B=s;}
	private:
	int _n;
	mat _X;
	Laplace *L;
	mat _B;
	mat _m;
	//vec _index;
	mat _s;
	RandomG::Random<boost::random::uniform_01<> > *runif;
	boost::math::normal_distribution<> *dn; 

};

